home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 12 / CU Amiga Magazine's Super CD-ROM 12 (1997)(EMAP Images)(GB)[!][issue 1997-07].iso / CUCD / Sound / MusicIn / psy.c < prev    next >
C/C++ Source or Header  |  1996-02-16  |  21KB  |  532 lines
  1. /**********************************************************************
  2. Copyright (c) 1991 MPEG/audio software simulation group, All Rights Reserved
  3. psy.c
  4. **********************************************************************/
  5. /**********************************************************************
  6.  * MPEG/audio coding/decoding software, work in progress              *
  7.  *   NOT for public distribution until verified and approved by the   *
  8.  *   MPEG/audio committee.  For further information, please contact   *
  9.  *   Davis Pan, 508-493-2241, e-mail: pan@3d.enet.dec.com             *
  10.  *                                                                    *
  11.  * VERSION 3.9                                                        *
  12.  *   changes made since last update:                                  *
  13.  *   date   programmers         comment                               *
  14.  * 2/25/91  Davis Pan           start of version 1.0 records          *
  15.  * 5/10/91  W. Joseph Carter    Ported to Macintosh and Unix.         *
  16.  * 7/10/91  Earle Jennings      Ported to MsDos.                      *
  17.  *                              replace of floats with FLOAT          *
  18.  * 2/11/92  W. Joseph Carter    Fixed mem_alloc() arg for "absthr".   *
  19.  * 8/07/92  Mike Coleman        Bug fix, read_absthr()                *
  20.   * 7/24/92  M. Iwadare          HANN window coefficients modified.    *
  21.  * 7/27/92  Masahiro Iwadare    Bug fix, FFT modification for Layer 3 *
  22.  * 7/27/92  Masahiro Iwadare    Bug fix, "new", "old", and "oldest"   *
  23.  *                              updates                               *
  24.  * 8/07/92  Mike Coleman        Bug fix, read_absthr()                *
  25.  * 14/07/95 Stephane Tavenard   automatic variables -> static (alloc) *
  26.  **********************************************************************/
  27.  
  28. #define ST_OPTIMIZE
  29.  
  30. #include "common.h"
  31. #include "encoder.h"
  32.  
  33. void psycho_anal(buffer,savebuf,chn,lay,snr32,sfreq)
  34. short int *buffer;
  35. short int savebuf[1056];
  36. int   chn, lay;
  37. FLOAT snr32[32];
  38. double sfreq;        /* to match prototype : float args are always double */
  39. {
  40.  unsigned int   i, j, k;
  41.  FLOAT          r_prime, phi_prime;
  42.  FLOAT          freq_mult, bval_lo, minthres, sum_energy;
  43.  double         tb, temp1, temp2, temp3;
  44.  
  45. /* The static variables "r", "phi_sav", "new", "old" and "oldest" have    */
  46. /* to be remembered for the unpredictability measure.  For "r" and        */
  47. /* "phi_sav", the first index from the left is the channel select and     */
  48. /* the second index is the "age" of the data.                             */
  49.  
  50.  static int     new = 0, old = 1, oldest = 0;
  51.  static int     init = 0, flush, sync_flush, syncsize, sfreq_idx;
  52.  
  53. /* The following static variables are constants.                           */
  54.  
  55.  static double  nmt = 5.5;
  56.  
  57.  static FLOAT   crit_band[27] = {0,  100,  200, 300, 400, 510, 630,  770,
  58.                                920, 1080, 1270,1480,1720,2000,2320, 2700,
  59.                               3150, 3700, 4400,5300,6400,7700,9500,12000,
  60.                              15500,25000,30000};
  61.  
  62.  static FLOAT   bmax[27] = {20.0, 20.0, 20.0, 20.0, 20.0, 17.0, 15.0,
  63.                             10.0,  7.0,  4.4,  4.5,  4.5,  4.5,  4.5,
  64.                              4.5,  4.5,  4.5,  4.5,  4.5,  4.5,  4.5,
  65.                              4.5,  4.5,  4.5,  3.5,  3.5,  3.5};
  66.  
  67. /* The following pointer variables point to large areas of memory         */
  68. /* dynamically allocated by the mem_alloc() function.  Dynamic memory     */
  69. /* allocation is used in order to avoid stack frame or data area          */
  70. /* overflow errors that otherwise would have occurred at compile time     */
  71. /* on the Macintosh computer.                                             */
  72.  
  73.  static FLOAT          *grouped_c, *grouped_e, *nb, *cb, *ecb, *bc;
  74.  static FLOAT          *wsamp_r, *wsamp_i, *phi, *energy;
  75.  static FLOAT          *c, *fthr;
  76.  static F32            *snrtmp;
  77.  
  78.  static int     *numlines;
  79.  static int     *partition;
  80.  static FLOAT   *cbval, *rnorm;
  81.  static FLOAT   *window;
  82.  static FLOAT   *absthr;
  83.  static double  *tmn;
  84.  static FCB     *s;
  85.  static FHBLK   *lthr;
  86.  static F2HBLK  *r, *phi_sav;
  87.  
  88. #ifdef ST_OPTIMIZE     /* ST 14/07/95 */
  89.  FLOAT  *sp;
  90.  FLOAT  *ecbj;
  91.  FLOAT  *cbj;
  92.  FLOAT  *grc;
  93.  FLOAT  *gre;
  94.  
  95.  if( !init ) {
  96. #endif
  97. /* These dynamic memory allocations simulate "automatic" variables        */
  98. /* placed on the stack.  For each mem_alloc() call here, there must be    */
  99. /* a corresponding mem_free() call at the end of this function.           */
  100.  
  101.  grouped_c = (FLOAT *) mem_alloc(sizeof(FCB), "grouped_c");
  102.  grouped_e = (FLOAT *) mem_alloc(sizeof(FCB), "grouped_e");
  103.  nb = (FLOAT *) mem_alloc(sizeof(FCB), "nb");
  104.  cb = (FLOAT *) mem_alloc(sizeof(FCB), "cb");
  105.  ecb = (FLOAT *) mem_alloc(sizeof(FCB), "ecb");
  106.  bc = (FLOAT *) mem_alloc(sizeof(FCB), "bc");
  107.  wsamp_r = (FLOAT *) mem_alloc(sizeof(FBLK), "wsamp_r");
  108.  wsamp_i = (FLOAT *) mem_alloc(sizeof(FBLK), "wsamp_i");
  109.  phi = (FLOAT *) mem_alloc(sizeof(FBLK), "phi");
  110.  energy = (FLOAT *) mem_alloc(sizeof(FBLK), "energy");
  111.  c = (FLOAT *) mem_alloc(sizeof(FHBLK), "c");
  112.  fthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "fthr");
  113.  snrtmp = (F32 *) mem_alloc(sizeof(F2_32), "snrtmp");
  114. #ifndef ST_OPTIMIZE      /* ST 14/07/95 */
  115.  if(init==0){
  116. #endif
  117.  
  118. /* These dynamic memory allocations simulate "static" variables placed    */
  119. /* in the data space.  Each mem_alloc() call here occurs only once at     */
  120. /* initialization time.  The mem_free() function must not be called.      */
  121.  
  122.      numlines = (int *) mem_alloc(sizeof(ICB), "numlines");
  123.      partition = (int *) mem_alloc(sizeof(IHBLK), "partition");
  124.      cbval = (FLOAT *) mem_alloc(sizeof(FCB), "cbval");
  125.      rnorm = (FLOAT *) mem_alloc(sizeof(FCB), "rnorm");
  126.      window = (FLOAT *) mem_alloc(sizeof(FBLK), "window");
  127.      absthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "absthr");
  128.      tmn = (double *) mem_alloc(sizeof(DCB), "tmn");
  129.      s = (FCB *) mem_alloc(sizeof(FCBCB), "s");
  130.      lthr = (FHBLK *) mem_alloc(sizeof(F2HBLK), "lthr");
  131.      r = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "r");
  132.      phi_sav = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "phi_sav");
  133.  
  134.      i = sfreq + 0.5;
  135.      switch(i){
  136.         case 32000: sfreq_idx = 0; break;
  137.         case 44100: sfreq_idx = 1; break;
  138.         case 48000: sfreq_idx = 2; break;
  139.         default:    printf("error, invalid sampling frequency: %d Hz\n",i);
  140.         exit(-1);
  141.      }
  142. /*     printf("absthr[][] sampling frequency index: %d\n",sfreq_idx); */
  143.      read_absthr(absthr, sfreq_idx);
  144.      if(lay==1){
  145.         flush = 384;
  146.         syncsize = 1024;
  147.         sync_flush = 576;
  148.      }
  149.      else {
  150.         flush = 384*3.0/2.0;
  151.         syncsize = 1056;
  152.         sync_flush = syncsize - flush;
  153.      }
  154. /* calculate HANN window coefficients */
  155. /*   for(i=0;i<BLKSIZE;i++)window[i]=0.5*(1-cos(2.0*PI*i/(BLKSIZE-1.0))); */
  156.      for(i=0;i<BLKSIZE;i++)window[i]=0.5*(1-cos(2.0*PI*(i-0.5)/BLKSIZE));
  157. /* reset states used in unpredictability measure */
  158.      for(i=0;i<HBLKSIZE;i++){
  159.         r[0][0][i]=r[1][0][i]=r[0][1][i]=r[1][1][i]=0;
  160.         phi_sav[0][0][i]=phi_sav[1][0][i]=0;
  161.         phi_sav[0][1][i]=phi_sav[1][1][i]=0;
  162.         lthr[0][i] = 60802371420160.0;
  163.         lthr[1][i] = 60802371420160.0;
  164.      }
  165. /*****************************************************************************
  166.  * Initialization: Compute the following constants for use later             *
  167.  *    partition[HBLKSIZE] = the partition number associated with each        *
  168.  *                          frequency line                                   *
  169.  *    cbval[CBANDS]       = the center (average) bark value of each          *
  170.  *                          partition                                        *
  171.  *    numlines[CBANDS]    = the number of frequency lines in each partition  *
  172.  *    tmn[CBANDS]         = tone masking noise                               *
  173.  *****************************************************************************/
  174. /* compute fft frequency multiplicand */
  175.      freq_mult = sfreq/BLKSIZE;
  176.  
  177. /* calculate fft frequency, then bval of each line (use fthr[] as tmp storage)*/
  178.      for(i=0;i<HBLKSIZE;i++){
  179.         temp1 = i*freq_mult;
  180.         j = 1;
  181.         while(temp1>crit_band[j])j++;
  182.         fthr[i]=j-1+(temp1-crit_band[j-1])/(crit_band[j]-crit_band[j-1]);
  183.      }
  184.      partition[0] = 0;
  185. /* temp2 is the counter of the number of frequency lines in each partition */
  186.      temp2 = 1;
  187.      cbval[0]=fthr[0];
  188.      bval_lo=fthr[0];
  189.      for(i=1;i<HBLKSIZE;i++){
  190.         if((fthr[i]-bval_lo)>0.33){
  191.            partition[i]=partition[i-1]+1;
  192.            cbval[partition[i-1]] = cbval[partition[i-1]]/temp2;
  193.            cbval[partition[i]] = fthr[i];
  194.            bval_lo = fthr[i];
  195.            numlines[partition[i-1]] = temp2;
  196.            temp2 = 1;
  197.         }
  198.         else {
  199.            partition[i]=partition[i-1];
  200.            cbval[partition[i]] += fthr[i];
  201.            temp2++;
  202.         }
  203.      }
  204.      numlines[partition[i-1]] = temp2;
  205.      cbval[partition[i-1]] = cbval[partition[i-1]]/temp2;
  206.  
  207. /************************************************************************
  208.  * Now compute the spreading function, s[j][i], the value of the spread-*
  209.  * ing function, centered at band j, for band i, store for later use    *
  210.  ************************************************************************/
  211.      for(j=0;j<CBANDS;j++){
  212.         for(i=0;i<CBANDS;i++){
  213.            temp1 = (cbval[i] - cbval[j])*1.05;
  214.            if(temp1>=0.5 && temp1<=2.5){
  215.               temp2 = temp1 - 0.5;
  216.               temp2 = 8.0 * (temp2*temp2 - 2.0 * temp2);
  217.            }
  218.            else temp2 = 0;
  219.            temp1 += 0.474;
  220.            temp3 = 15.811389+7.5*temp1-17.5*sqrt((double) (1.0+temp1*temp1));
  221.            if(temp3 <= -100) s[i][j] = 0;
  222.            else {
  223.               temp3 = (temp2 + temp3)*LN_TO_LOG10;
  224.               s[i][j] = exp(temp3);
  225.            }
  226.         }
  227.      }
  228.  
  229.   /* Calculate Tone Masking Noise values */
  230.      for(j=0;j<CBANDS;j++){
  231.         temp1 = 15.5 + cbval[j];
  232.         tmn[j] = (temp1>24.5) ? temp1 : 24.5;
  233.   /* Calculate normalization factors for the net spreading functions */
  234.         rnorm[j] = 0;
  235.         for(i=0;i<CBANDS;i++){
  236.            rnorm[j] += s[j][i];
  237.         }
  238.      }
  239.      init++;
  240.  }
  241.  
  242. /************************* End of Initialization *****************************/
  243.  switch(lay) {
  244.   case 1:
  245.   case 2:
  246.      for(i=0; i<lay; i++){
  247. /*****************************************************************************
  248.  * Net offset is 480 samples (1056-576) for layer 2; this is because one must*
  249.  * stagger input data by 256 samples to synchronize psychoacoustic model with*
  250.  * filter bank outputs, then stagger so that center of 1024 FFT window lines *
  251.  * up with center of 576 "new" audio samples.                                *
  252.  *                                                                           *
  253.  * For layer 1, the input data still needs to be staggered by 256 samples,   *
  254.  * then it must be staggered again so that the 384 "new" samples are centered*
  255.  * in the 1024 FFT window.  The net offset is then 576 and you need 448 "new"*
  256.  * samples for each iteration to keep the 384 samples of interest centered   *
  257.  *****************************************************************************/
  258.         for(j=0; j<syncsize; j++){
  259.            if(j<(sync_flush))savebuf[j] = savebuf[j+flush];
  260.            else savebuf[j] = *buffer++;
  261.            if(j<BLKSIZE){
  262. /**window data with HANN window***********************************************/
  263.               wsamp_r[j] = window[j]*((FLOAT) savebuf[j]);
  264.               wsamp_i[j] = 0;
  265.            }
  266.         }
  267. /**Compute FFT****************************************************************/
  268.         fft(wsamp_r,wsamp_i,energy,phi,1024);
  269. /*****************************************************************************
  270.  * calculate the unpredictability measure, given energy[f] and phi[f]        *
  271.  *****************************************************************************/
  272. /*only update data "age" pointers after you are done with both channels      */
  273. /*for layer 1 computations, for the layer 2 double computations, the pointers*/
  274. /*are reset automatically on the second pass                                 */
  275.          if(lay==2 || (lay==1 && chn==0) ){
  276.            if(new==0){new = 1; oldest = 1;}
  277.            else {new = 0; oldest = 0;}
  278.            if(old==0)old = 1; else old = 0;
  279.         }
  280. #ifdef ST_OPTIMIZE /* ST 14/08/95 */
  281.         {
  282.            FLOAT *r_old = &r[chn][old][0];
  283.            FLOAT *r_oldest = &r[chn][oldest][0];
  284.            FLOAT *r_new = &r[chn][new][0];
  285.            FLOAT *phi_old = &phi_sav[chn][old][0];
  286.            FLOAT *phi_oldest = &phi_sav[chn][oldest][0];
  287.            FLOAT *phi_new = &phi_sav[chn][new][0];
  288.            FLOAT *en = &energy[0];
  289.            FLOAT *ph = &phi[0];
  290.            FLOAT *cj = &c[0];
  291.  
  292.            for( j=0; j<HBLKSIZE; j++ ) {
  293.               r_prime = 2 * *r_old++ - *r_oldest++;
  294.               phi_prime = 2 * *phi_old++ - *phi_oldest++;
  295.               *r_new = sqrt( (double)*en++ );
  296.               *phi_new = *ph++;
  297.               temp1 = *r_new * cos((double)*phi_new) - r_prime * cos((double)phi_prime);
  298.               temp2 = *r_new * sin((double)*phi_new++) - r_prime * sin((double)phi_prime);
  299.               temp3 = *r_new++ + fabs((double)r_prime);
  300.               if( temp3 != 0 ) *cj++ = sqrt( temp1*temp1 + temp2*temp2 ) / temp3;
  301.               else *cj++ = 0;
  302.            }
  303.         }
  304. #else
  305.         for(j=0; j<HBLKSIZE; j++){
  306.            r_prime = 2.0 * r[chn][old][j] - r[chn][oldest][j];
  307.            phi_prime = 2.0 * phi_sav[chn][old][j] - phi_sav[chn][oldest][j];
  308.            r[chn][new][j] = sqrt((double) energy[j]);
  309.            phi_sav[chn][new][j] = phi[j];
  310. temp1=r[chn][new][j] * cos((double) phi[j]) - r_prime * cos((double) phi_prime);
  311. temp2=r[chn][new][j] * sin((double) phi[j]) - r_prime * sin((double) phi_prime);
  312.            temp3=r[chn][new][j] + fabs((double)r_prime);
  313.            if(temp3 != 0)c[j]=sqrt(temp1*temp1+temp2*temp2)/temp3;
  314.            else c[j] = 0;
  315.         }
  316. #endif
  317. /*****************************************************************************
  318.  * Calculate the grouped, energy-weighted, unpredictability measure,         *
  319.  * grouped_c[], and the grouped energy. grouped_e[]                          *
  320.  *****************************************************************************/
  321. #ifdef ST_OPTIMIZE /* ST 14/07/95 */
  322.         memset( grouped_c, 0, sizeof(FCB) );
  323.         memset( grouped_e, 0, sizeof(FCB) );
  324. #else
  325.         for(j=1;j<CBANDS;j++){
  326.            grouped_e[j] = 0;
  327.            grouped_c[j] = 0;
  328.         }
  329. #endif
  330.         grouped_e[0] = energy[0];
  331.         grouped_c[0] = energy[0]*c[0];
  332.         for(j=1;j<HBLKSIZE;j++){
  333.            grouped_e[partition[j]] += energy[j];
  334.            grouped_c[partition[j]] += energy[j]*c[j];
  335.         }
  336. /*****************************************************************************
  337.  * convolve the grouped energy-weighted unpredictability measure             *
  338.  * and the grouped energy with the spreading function, s[j][k]               *
  339.  *****************************************************************************/
  340. #ifdef ST_OPTIMIZE /* ST 14/07/95 */
  341.         sp = &s[0][0];
  342.         ecbj = ecb;
  343.         cbj = cb;
  344.         for( j=0; j<CBANDS; j++ ) {
  345.            *ecbj = 0;
  346.            *cbj = 0;
  347.            gre = grouped_e;
  348.            grc = grouped_c;
  349.            for( k=0; k<CBANDS; k++ ) {
  350.               if( *sp != 0.0 ) {
  351.                  *ecbj += *sp * *gre;
  352.                  *cbj += *sp * *grc;
  353.               }
  354.               gre++;
  355.               grc++;
  356.               sp++;
  357.            }
  358.            if( *ecbj != 0 ) {
  359.               *cbj = *cbj / *ecbj;
  360.               if( *cbj > 0.5 ) {
  361.                  *cbj = 0.5;
  362.                  tb = 0.0;
  363.               }
  364.               else tb = -0.434294482*log((double)*cbj)-0.301029996;
  365.            }
  366.            else {
  367.               *cbj = 0.05;
  368.               tb  = 1.0;
  369.            }
  370.            ecbj++;
  371.            cbj++;
  372.  
  373.            bc[j] = tmn[j]*tb + nmt*(1.0-tb);
  374.            k = cbval[j] + 0.5;
  375.            bc[j] = (bc[j] > bmax[k]) ? bc[j] : bmax[k];
  376.            bc[j] = exp((double) -bc[j]*LN_TO_LOG10);
  377.         }
  378. #else
  379.         for(j=0;j<CBANDS;j++){
  380.            ecb[j] = 0;
  381.            cb[j] = 0;
  382.            for(k=0;k<CBANDS;k++){
  383.               if(s[j][k] != 0.0){
  384.                  ecb[j] += s[j][k]*grouped_e[k];
  385.                  cb[j] += s[j][k]*grouped_c[k];
  386.               }
  387.            }
  388.            if(ecb[j] !=0)cb[j] = cb[j]/ecb[j];
  389.            else cb[j] = 0;
  390.         }
  391. /*****************************************************************************
  392.  * Calculate the required SNR for each of the frequency partitions           *
  393.  *         this whole section can be accomplished by a table lookup          *
  394.  *****************************************************************************/
  395.         for(j=0;j<CBANDS;j++){
  396.            if(cb[j]<.05)cb[j]=0.05;
  397.            else if(cb[j]>.5)cb[j]=0.5;
  398.            tb = -0.434294482*log((double) cb[j])-0.301029996;
  399.            bc[j] = tmn[j]*tb + nmt*(1.0-tb);
  400.            k = cbval[j] + 0.5;
  401.            bc[j] = (bc[j] > bmax[k]) ? bc[j] : bmax[k];
  402.            bc[j] = exp((double) -bc[j]*LN_TO_LOG10);
  403.         }
  404. #endif
  405. /*****************************************************************************
  406.  * Calculate the permissible noise energy level in each of the frequency     *
  407.  * partitions. Include absolute threshold and pre-echo controls              *
  408.  *         this whole section can be accomplished by a table lookup          *
  409.  *****************************************************************************/
  410.         for(j=0;j<CBANDS;j++)
  411.            if(rnorm[j] && numlines[j])
  412.               nb[j] = ecb[j]*bc[j]/(rnorm[j]*numlines[j]);
  413.            else nb[j] = 0;
  414.         for(j=0;j<HBLKSIZE;j++){
  415. /*temp1 is the preliminary threshold */
  416.            temp1=nb[partition[j]];
  417.            temp1=(temp1>absthr[j])?temp1:absthr[j];
  418. /*do not use pre-echo control for layer 2 because it may do bad things to the*/
  419. /*  MUSICAM bit allocation algorithm                                         */
  420.            if(lay==1){
  421.               fthr[j] = (temp1 < lthr[chn][j]) ? temp1 : lthr[chn][j];
  422.               temp2 = temp1 * 0.00316;
  423.               fthr[j] = (temp2 > fthr[j]) ? temp2 : fthr[j];
  424.            }
  425.            else fthr[j] = temp1;
  426.            lthr[chn][j] = LXMIN*temp1;
  427.         }
  428. /*****************************************************************************
  429.  * Translate the 512 threshold values to the 32 filter bands of the coder    *
  430.  *****************************************************************************/
  431.         for(j=0;j<193;j += 16){
  432.            minthres = 60802371420160.0;
  433.            sum_energy = 0.0;
  434.            for(k=0;k<17;k++){
  435.               if(minthres>fthr[j+k])minthres = fthr[j+k];
  436.               sum_energy += energy[j+k];
  437.            }
  438.            snrtmp[i][j/16] = sum_energy/(minthres * 17.0);
  439.            snrtmp[i][j/16] = 4.342944819 * log((double)snrtmp[i][j/16]);
  440.         }
  441.         for(j=208;j<(HBLKSIZE-1);j += 16){
  442.            minthres = 0.0;
  443.            sum_energy = 0.0;
  444.            for(k=0;k<17;k++){
  445.               minthres += fthr[j+k];
  446.               sum_energy += energy[j+k];
  447.            }
  448.            snrtmp[i][j/16] = sum_energy/minthres;
  449.            snrtmp[i][j/16] = 4.342944819 * log((double)snrtmp[i][j/16]);
  450.         }
  451. /*****************************************************************************
  452.  * End of Psychoacuostic calculation loop                                    *
  453.  *****************************************************************************/
  454.      }
  455.      for(i=0; i<32; i++){
  456.         if(lay==2)
  457.            snr32[i]=(snrtmp[0][i]>snrtmp[1][i])?snrtmp[0][i]:snrtmp[1][i];
  458.         else snr32[i]=snrtmp[0][i];
  459.      }
  460.      break;
  461.   case 3:
  462.      printf("layer 3 is not currently supported\n");
  463.      break;
  464.   default:
  465.      printf("error, invalid MPEG/audio coding layer: %d\n",lay);
  466.  }
  467.  
  468. #ifndef ST_OPTIMIZE      /* ST 14/07/95 */
  469. /* These mem_free() calls must correspond with the mem_alloc() calls     */
  470. /* used at the beginning of this function to simulate "automatic"        */
  471. /* variables placed on the stack.                                        */
  472.  
  473.  mem_free((void **) &grouped_c);
  474.  mem_free((void **) &grouped_e);
  475.  mem_free((void **) &nb);
  476.  mem_free((void **) &cb);
  477.  mem_free((void **) &ecb);
  478.  mem_free((void **) &bc);
  479.  mem_free((void **) &wsamp_r);
  480.  mem_free((void **) &wsamp_i);
  481.  mem_free((void **) &phi);
  482.  mem_free((void **) &energy);
  483.  mem_free((void **) &c);
  484.  mem_free((void **) &fthr);
  485.  mem_free((void **) &snrtmp);
  486. #endif
  487.  
  488. }
  489.  
  490. /******************************************************************************
  491. routine to read in absthr table from a file.
  492. ******************************************************************************/
  493.  
  494. void read_absthr(absthr, table)
  495. FLOAT *absthr;
  496. int table;
  497. {
  498.  FILE *fp;
  499.  long j,index;
  500.  float a;
  501.  char t[80];
  502.  char ta[16];
  503.  
  504.  strcpy( ta, "absthr_0" );
  505.  
  506.  switch(table){
  507.     case 0 : ta[7] = '0';
  508.              break;
  509.     case 1 : ta[7] = '1';
  510.              break;
  511.     case 2 : ta[7] = '2';
  512.              break;
  513.     default : printf("absthr table: Not valid table number\n");
  514.  }
  515.  if(!(fp = OpenTableFile(ta) ) ){
  516.     printf("Please check %s table\n", ta);
  517.     exit(1);
  518.  }
  519.  fgets(t, 150, fp);
  520.  sscanf(t, "table %ld", &index);
  521.  if(index != table){
  522.     printf("error in absthr table %s",ta);
  523.     exit(1);
  524.  }
  525.  for(j=0; j<HBLKSIZE; j++){
  526.     fgets(t,80,fp);
  527.     sscanf(t,"%f", &a);
  528.     absthr[j] =  a;
  529.  }
  530.  fclose(fp);
  531. }
  532.